function[stimes, sval] = stairsum(jmptimes, fval)
%
% STAIRSUM Add piecewise constant (stair) functions.  
%   The functions are stored in two matrices containing the jumps points and
%   function values respectively. If they have different number of jumps,
%   each vector is padded with the last value to have the same length.
%   Functions are assumed to be right-continuous. The superposition
%   process is also a piecewise constant function.
%
% [stimes, sval] = stairsum(jmptimes, fval)
%
% Inputs:
%        jmptimes - a matrix of the jump times stored columnwise
%        fval - a matrix of the function values at the jump points 
%          stored columnwise 
%
% Outputs:
%        stimes - a column vector of the jump times of the
%          superposition process 
%        sval - a column vector of values of the superposition
%          process at the jump points
%
% See also STAIRCUT, STAIRINTEGR.

% Authors: R.Gaigalas, I.Kaj
% v2.0 18-Oct-05

  % starting with R13 we may use numel here
  ntotal = size(jmptimes, 1)*size(jmptimes, 2);
  
  % merge all jump times into one column vector
  all_times = reshape(jmptimes, ntotal, 1);
  
  % "differentiate": find increments of the processes and 
  % merge into one column vector 
  all_jumps = reshape([fval(1, :); diff(fval)], ntotal, 1);

  % sort the "delta" functions at the jump times
  ppsum = sortrows([all_times all_jumps], 1);
  stimes = ppsum(:, 1);
  sval = cumsum(ppsum(:, 2));
  
